#install.packages('devtools')
#devtools::install_github("jinworks/CellChat")
library(CellChat)
library(tidyverse)
library(ggalluvial)
# install.packages("anndata")
library(anndata)
library(anndata)
library(Seurat)
library(patchwork)
library(uwot)
library(reticulate)
library(glmGamPoi)
# BiocManager::install("glmGamPoi")
library(umap)
# install.packages("umap")
library(ComplexHeatmap)
#install.packages('NMF')
library(NMF)
library(dplyr)
# BiocManager::install("SeuratData")
library(SeuratData)
library(ggplot2)
library(svglite)
# install.packages('wordcloud','wordcloud2')
library(wordcloud)
library(wordcloud2)
# BiocManager::install("tm")
library(tm)
#devtools::install_github("jokergoo/circlize")
library(circlize)
#devtools::install_github("jokergoo/ComplexHeatmap")

#source('functional.R')#千萬不要加上
# setwd("~/cellchat")
# rm(list=ls()) #清空所有变量
options(stringsAsFactors = FALSE)#输入数据不自动转换成因子（防止数据格式错误）

#(A) Starting from a count data matrix
# load("data_humanSkin_CellChat.rda")
# data.input = data_humanSkin$data
# meta = data_humanSkin$meta
# cell.use = rownames(meta)[meta$condition == "LS"] # extract the c
# data.input = data.input[, cell.use]
# meta = meta[cell.use, ]


#(B) Starting from a Seurat object
#data.input <- GetAssayData(seurat_object, assay = "RNA", slot = "data") # normalized data matrix
#labels <- Idents(seurat_object)
#meta <- data.frame(group = labels, row.names = names(labels)) # create a dataframe of the cell labels


#⭐合并分析####
#1. 读入数据####
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")
#T_cell
obj1 <- readRDS("./data/temp/T_cluster_id_test_1.rds")
levels(x = obj1)
obj1 <- subset(obj1,idents= c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm") )

#Duct_epithelial_cell
obj2 <- readRDS("./data/temp/Duct_epithelial_cell_cluster_id_test.rds")
levels(x = obj2)

#合并
obj_all <- merge(obj1,obj2)
#labels细胞分群名，其他分群也这样修改⭐
obj_all$labels = as.character(Idents(obj_all))
levels(x = obj_all)
table(obj_all@meta.data[["labels"]])
rm(list = c("obj1","obj2","all"))

obj_all <- saveRDS("./data/temp/T-Duct_cellchat.all_obj_all.rds")

#2. 转化为cellchat数据格式####
#(A)从矩阵创建
#取数据all@assay$RNA@data
# data.input <- GetAssayData(obj_all, assay = "RNA", slot = "data") # normalized data matrix
# meta <- data.frame(group = labels, row.names = names(labels)) # create a dataframe of the cell labels
# # meta <- data.frame(group = labels, row.names = names(active.idents)) # create a dataframe of the cell labels
# cellchat <- createCellChat(object = data.input, meta = meta, group.by = "RNA")#group.by = "labels"

#(B)从Seurat创建
cellchat <- createCellChat(object = obj_all, group.by = "labels", assay = "RNA")


#3. 设置参考数据库####
# 根据分析数据的物种，可选CellChatDB.human, 或者 CellChatDB.mouse 。通过showDatabaseCategory函数可以查看该数据库的情况
CellChatDB <- CellChatDB.human
showDatabaseCategory(CellChatDB)
# 展示参考数据库
dplyr::glimpse(CellChatDB$interaction)
# (A)使用默认所有的参考数据库
CellChatDB.use <- CellChatDB # 

# (B)使用CellChatDB的子集
# CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")# use Secreted Signaling

# 合并保存到cellchat内
cellchat@DB <- CellChatDB.use


# 对信号基因的表达数据取子集以节省计算成本
cellchat <- subsetData(cellchat)                         # 即使使用整个数据库，此步骤也是必要的
future::plan("multisession", workers = 10)                # 运行线程数
cellchat <- identifyOverExpressedGenes(cellchat)         #识别过表达基因
cellchat <- identifyOverExpressedInteractions(cellchat)  #识别过表达配体受体对

#细胞通信网络的推理
cellchat <- computeCommunProb(cellchat, type = "triMean") #triMean用于计算每个细胞组的平均基因表达。


ptm <-  Sys.time() 
execution.time = Sys.time() - ptm
print(as.numeric(execution.time, units = "secs"))

# cellchat <- filterCommunication(cellchat, min.cells = 10)
# cellchat <- computeCommunProbPathway(cellchat)

#4. 推断细胞通讯网络####
#使用表达值推测细胞互作的概率，该步骤相对较耗时一些。
cellchat <- computeCommunProb(cellchat, raw.use = TRUE, population.size = TRUE) 
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
# 注1：raw.use = TRUE 表示使用raw数据，而不使用上一步projectData后的结果。
# 注2：在假设细胞数较多的群 往往比 细胞数较少的群发送更强的信号的前提下，当population.size = TRUE时候，CellChat可以在概率计算中考虑每个细胞群中细胞比例的影响。
cellchat <- filterCommunication(cellchat, min.cells = 10)



#5. 提取 保存结果####
#all the inferred cell-cell communications at the level of ligands/receptors
df.net <- subsetCommunication(cellchat)
write.csv(df.net, "./data/output/T-Duct_cellchat.all.csv")

#获取显著的结果
df.net1 <- subsetCommunication(cellchat,slot.name = "netP")

#获取起始和结束的关系
levels(cellchat@idents)
df.net2 <- subsetCommunication(cellchat, sources.use = c("Epi"), targets.use = c("Fibroblast" ,"T")) 

#获取信号通路
df.net3 <- subsetCommunication(cellchat, signaling = c("CCL", "TGFb"))

#6. 计算cell-cell communication####
#计算每个信号通路相关的所有配体-受体相互作用的通信结果，结存存放在net 和 netP中 。
cellchat <- computeCommunProbPathway(cellchat)

#计算整合的细胞类型之间通信结果
cellchat <- aggregateNet(cellchat)

#耗时
execution.time = Sys.time() - ptm
print(as.numeric(execution.time, units = "secs"))

#保存文件
saveRDS(cellchat, file = "./data/temp/T-Duct_cellchat.all.rds")
#20240612@@@@


#7. 可视化####
cellchat <- readRDS(file = "./data/temp/T-Duct_cellchat.all.rds")
## Access all the signaling pathways showing significant communications
pathways.show.all <- cellchat@netP$pathways
pathways.show.all
# [1] "MHC-I"     "MIF"       "LAMININ"   "COLLAGEN"  "CLEC"      "SPP1"      "MK"       
# [8] "CD99"      "APP"       "FN1"       "GALECTIN"  "PARs"      "ADGRE5"    "CDH"      
# [15] "MHC-II"    "CDH1"      "CD46"      "LCK"       "SEMA4"     "VISFATIN"  "JAM"      
# [22] "EGF"       "NECTIN"    "TGFb"      "TIGIT"     "DESMOSOME" "NOTCH"     "CEACAM"   
# [29] "GRN"       "EPHA"      "AGRN"      "VTN"       "OCLN"      "THBS"      "HSPG"     
# [36] "GDF"       "ITGB2"     "TNF"       "ALCAM"     "CD6"       "CADM"      "CD96"     
# [43] "MPZ"       "SELPLG"    "VCAM"      "SEMA3"     "EPHB"      "CD226"     "INSULIN"  
# [50] "CXCL"      "PVR"       "IFN-II"    "ANGPTL"    "CCL"       "SEMA6"     "SELL"     
# [57] "TRAIL"     "BMP"       "ICAM"      "WNT"       "FLT3"      "SEMA7"     "TWEAK"  

##1. 单个信号通路####
# select one pathway
pathways.show <- c("TIGIT") #⭐T和Duct的交互
pathways.show <- c("EGF")#内部互相促进
pathways.show <- c("TNF")#⭐ OK
pathways.show <- c("SPP1")#
pathways.show <- c("LAMININ")#配对中有一个分别作用T和Duct，与tumor匹配
pathways.show <- c("MHC-II")#HLA-DMA - CD4,与tumor匹配


pathways.show <- c("COLLAGEN")
pathways.show <- c("CD99")#仅一个配对
#舍弃

pathways.show <- c("CDH")
pathways.show <- c("MHC-I") #
pathways.show <- c("ADGRE5")#仅一个配对
pathways.show <- c("PARs")#无特色
pathways.show <- c("CD99")#仅一个配对
pathways.show <- c("APP")#仅一个配对
pathways.show <- c("COLLAGEN")
pathways.show <- c("MIF")
pathways.show <- c("FN1")
pathways.show <- c("TGFb")
pathways.show <- c("NOTCH")#
pathways.show <- c("WNT")#
pathways.show <- c("LAMININ")#
pathways.show <- c("EPHA")


###1. 环状图####
par(mfrow=c(1,1))
p1 <- netVisual_aggregate(cellchat, signaling = pathways.show, 
                          layout = "circle", 
                          color.use = NULL, 
                          #label.edge= T,#显示贡献值
                          sources.use = NULL, 
                          targets.use = NULL, 
                          idents.use = NULL)

###2. 层次聚类图####
vertex.receiver = seq(1,9)
p2 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout ="hierarchy", vertex.receiver = vertex.receiver)
# 左半部分是自分泌相关信号，自己释放的信号作用于自己
# 右半部分就是展示的旁分泌信号
# 链接：https://www.jianshu.com/p/cf79a2b7703b

###3. 和弦图####
par(mfrow=c(1,1))
p3 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout ="chord")

par(mfrow=c(1,1))
group.cellType <- c(rep("CD4", 4), rep("CD8", 5),rep("Duct", 9), rep("NKT", 1) ) #grouping cell clusters into fibroblast, DC and TC cells
names(group.cellType) <- levels(cellchat@idents)
p4 <- netVisual_chord_cell(cellchat, signaling = pathways.show, group =group.cellType, title.name = paste0(pathways.show, " signaling network"))
p4
###4. 热图####
par(mfrow=c(1,1))
p5 <- netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
p5
###5. 柱状图#### 
# 配受体贡献列表
p6 <- netAnalysis_contribution(cellchat, signaling = pathways.show)
p6
配受体贡献列表配受体贡献列表
pairLR <- extractEnrichedLR(cellchat, signaling = pathways.show, geneLR.return = FALSE)
pairLR

###5. 配受体作用情况####
LR.show <- pairLR[1,]
p71 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[2,]
p72 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[3,]
p73 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[4,]
p74 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[5,]
p75 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[6,]
p76 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[7,]
p77 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[8,]
p78 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[9,]
p79 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[10,]
p710 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[11,]
p711 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[12,]
p712 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[13,]
p713 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")




pdf("./data/output/T-Duct_cellchat_MHC-II.pdf",width =10 ,height = 10)
p1
p2
p3
p4
p5
p6
p71
p72
p73
p74
p75
p76
p77
p78
p79
p710
p711
p712
p713
dev.off()


##2. 气泡图-出发结束####
# sources.use, targets.use 对应序号 
levels(cellchat@idents)
# [1] "CD4_Tem"  "CD4_Th"   "CD4_Tn"   "CD4_Treg" "CD8_Tc"   "CD8_Te"   "CD8_Tem"  "CD8_Tex"  "CD8_Trm" 
# [10] "Group_1"  "Group_2"  "Group_3"  "Group_4"  "Group_5"  "Group_6"  "Group_7"  "Group_8"  "Group_9" 
# [19] "NKT" 

###1. 所有的信号通路####
# T出发
p1 <- netVisual_bubble(cellchat, sources.use = c(1:9,19), targets.use = c(10:18), remove.isolate = FALSE)#去掉空白行和空白列
# Duct出发
p2 <- netVisual_bubble(cellchat, sources.use = c(10:18), targets.use = c(1:9,19), remove.isolate = FALSE)
pdf("./data/output/T-Duct_cellchat_All_气泡图.pdf",width =15 ,height = 10)
p1
p2
dev.off()

###2. 看指定信号通路####
pairLR.use <- extractEnrichedLR(cellchat, signaling = c("TIGIT","EGF","TNF","SPP1","MIF","MHC-II"))

netVisual_bubble(cellchat, sources.use = c(1:9,19), targets.use = c(10:18), pairLR.use = pairLR.use, remove.isolate = TRUE)
netVisual_bubble(cellchat, sources.use = c(10:18), targets.use = c(1:9,19),  
                 # pairLR.use = pairLR.use, 
                 remove.isolate = TRUE)

###3. 看内部信号通路####
#CD4内部 CD4Treg起点
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(1:4), 
                 #signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), 
                 remove.isolate = FALSE)
#CD4内部 CD4Treg终点
netVisual_bubble(cellchat, targets.use = 4, sources.use = c(1:4), 
                 # signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), 
                 remove.isolate = FALSE)
#CD8内部 CD8Tex起点
netVisual_bubble(cellchat, sources.use = 8, targets.use = c(6:10), 
                 signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), remove.isolate = FALSE)
#CD8内部 CD8Tex终点
netVisual_bubble(cellchat, targets.use = 8, sources.use = c(6:10), 
                 signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), remove.isolate = FALSE)

#列出所有与细胞相关的配体受体⭐
netVisual_bubble(cellchat, targets.use = c("CD4_Treg","CD8_Tex","NKT"), pairLR.use = pairLR.use, 
                 remove.isolate = TRUE, sort.by.target = T)
netVisual_bubble(cellchat, sources.use =c("CD4_Treg","CD8_Tex","NKT"), pairLR.use = pairLR.use, remove.isolate = TRUE,sort.by.source = T)

netVisual_bubble(cellchat, 
                 sources.use = c("CD4_Treg","CD8_Tex","NKT"), 
                 targets.use = c("CD4_Treg","CD8_Tex","NKT"), 
                 pairLR.use = pairLR.use, 
                 remove.isolate = TRUE, 
                 sort.by.source = T, 
                 sort.by.target = T, 
                 sort.by.source.priority = F)#排序优先级

##3. 和弦图####
#基因和信号通路关系展示
netVisual_chord_gene(cellchat, sources.use = 4, targets.use = c(1:4), lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = 8, targets.use = c(5:10), lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = 10, targets.use = c(5:10), lab.cex = 0.5,legend.pos.y = 30)

netVisual_chord_gene(cellchat, sources.use = c(1:4), targets.use = 10, lab.cex =0.5,legend.pos.y = 40)
netVisual_chord_gene(cellchat, sources.use = c(5:10), targets.use = 8, lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = c(5:10), targets.use = 10, lab.cex = 0.5,legend.pos.y = 30)

# show all the significant interactions (L-R pairs) associated with certain signaling pathways
# netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("CCL","CXCL"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c("CD4_Treg","CD8_Tex","NKT"), targets.use = c("CD4_Treg","CD8_Tex","NKT"), signaling =  c("TIGIT","EGF","TNF","SPP1","LAMININ","MHC-II"),legend.pos.x = 8, lab.cex = 0.5)

# show all the significant signaling pathways from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), slot.name = "netP", legend.pos.x = 10)

##4. 配受体小提琴图####
pathways.show <- c("TIGIT")

#通路配受体表达强度，小提琴图
plotGeneExpression(cellchat, signaling = pathways.show, enriched.only =TRUE)

#计算并可视化网络中心性得分
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")# the slot 'netP' means the inferred intercellular communication network of signaling pathways

#使用热图将计算的中心性得分可视化，识别细胞群的主要信号作用
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

#在2D空间中可视化主要源和目标
netAnalysis_signalingRole_scatter(cellchat, signaling = NULL)

#识别对某些细胞组的传出或传入信号贡献最大的信号
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")
ht1 + ht2

##5. 通信模式####
#(A)识别和可视化分泌细胞的outgoing通信模式
# 推断图案的数量。
selectK(cellchat, pattern ="outgoing")
#当传出模式的数量为3时，Cophenetic和Silhouette值都开始突然下降。
nPatterns = 6 #选择合适的分组数
cellchat <- identifyCommunicationPatterns(cellchat, pattern ="outgoing", k = nPatterns)
p1 <- netAnalysis_river(cellchat, pattern ="outgoing")# river plot
p2 <- netAnalysis_dot(cellchat, pattern ="outgoing")# dot plot

#(B)识别和可视化目标细胞的incoming通信模式
selectK(cellchat, pattern = "incoming")
nPatterns = 7 #
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns)
p3 <- netAnalysis_river(cellchat, pattern = "incoming")# river plot
p4 <- netAnalysis_dot(cellchat, pattern = "incoming")# dot plot

pdf("./data/output/T-Duct_cellchat_All_通信模式.pdf",width =8 ,height = 7)
p1
p2
p3
p4
dev.off()


##6. 相似性####
###1. 功能相似性Functional similarity####
cellchat <- computeNetSimilarity(cellchat, type = "functional")
# cellchat <- netEmbedding(cellchat, type = "functional")
#install through "pip install umap-learn"
#或者下面三條命令都可以運行
cellchat <- netEmbedding(cellchat, umap.method='uwot',type ="functional")
cellchat <- netEmbedding(cellchat, umap.method='umap-learn',type ="functional")
cellchat <- netClustering(cellchat, type = "functional",do.parallel=FALSE)#加上"do.parallel=FALSE"
p1 <- netVisual_embedding(cellchat, type = "functional", label.size =3.5)
p2 <- netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)

###2. 结构相似性Structure similarity####
cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")#或者下面這條命令
cellchat <- netEmbedding(cellchat, umap.method='uwot',type = "structural")
cellchat <- netClustering(cellchat, type = "structural",do.parallel=FALSE)#加上"do.parallel=FALSE"
p3 <- netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
p4 <- netVisual_embeddingZoomIn(cellchat, type = "structural", nCol = 2)

pdf("./data/output/T-Duct_cellchat_All_结构功能相似性.pdf",width =4 ,height = 3.5)
p1
p2
p3
p4
dev.off()

# saveRDS(cellchat, file = "cellchat_humanSkin_LS.rds")


#⭐ tumor####
##1. 读入数据####
obj_all <- saveRDS("./data/temp/T-Duct_cellchat.all_obj_all.rds")
obj_tumor <- subset(obj_all, tech=="Tumor")

##2. 转化为cellchat数据格式####
cellchat <- createCellChat(object = obj_tumor, group.by = "labels", assay = "RNA")
cellchat <- createCellChat(object = obj_normal, group.by = "labels", assay = "RNA")

##3. 设置参考数据库####
# 根据分析数据的物种，可选CellChatDB.human, 或者 CellChatDB.mouse 。通过showDatabaseCategory函数可以查看该数据库的情况
CellChatDB <- CellChatDB.human
showDatabaseCategory(CellChatDB)
# 展示参考数据库
dplyr::glimpse(CellChatDB$interaction)
# 使用默认所有的参考数据库
CellChatDB.use <- CellChatDB # 
# 合并保存到cellchat内
cellchat@DB <- CellChatDB.use

# 对信号基因的表达数据取子集以节省计算成本
cellchat <- subsetData(cellchat)                         # 即使使用整个数据库，此步骤也是必要的
future::plan("multisession", workers = 10)                # 运行线程数
cellchat <- identifyOverExpressedGenes(cellchat)         #识别过表达基因
cellchat <- identifyOverExpressedInteractions(cellchat)  #识别过表达配体受体对

#细胞通信网络的推理
cellchat <- computeCommunProb(cellchat, type = "triMean") #triMean用于计算每个细胞组的平均基因表达。

# cellchat <- filterCommunication(cellchat, min.cells = 10)
# cellchat <- computeCommunProbPathway(cellchat)

##4. 推断细胞通讯网络####
#使用表达值推测细胞互作的概率，该步骤相对较耗时一些。
cellchat <- computeCommunProb(cellchat, raw.use = TRUE, population.size = TRUE) 
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
# 注1：raw.use = TRUE 表示使用raw数据，而不使用上一步projectData后的结果。
# 注2：在假设细胞数较多的群 往往比 细胞数较少的群发送更强的信号的前提下，当population.size = TRUE时候，CellChat可以在概率计算中考虑每个细胞群中细胞比例的影响。
cellchat <- filterCommunication(cellchat, min.cells = 10)



##5. 提取 保存结果####
#all the inferred cell-cell communications at the level of ligands/receptors
df.net <- subsetCommunication(cellchat)
write.csv(df.net, "./data/output/T-Duct_cellchat.all_tumor.csv")

#获取显著的结果
df.net1 <- subsetCommunication(cellchat,slot.name = "netP")

#获取起始和结束的关系
levels(cellchat@idents)
df.net2 <- subsetCommunication(cellchat, sources.use = c("Epi"), targets.use = c("Fibroblast" ,"T")) 

#获取信号通路
df.net3 <- subsetCommunication(cellchat, signaling = c("CCL", "TGFb"))

##6. 计算cell-cell communication####
#计算每个信号通路相关的所有配体-受体相互作用的通信结果，结存存放在net 和 netP中 。
cellchat <- computeCommunProbPathway(cellchat)

#计算整合的细胞类型之间通信结果
cellchat <- aggregateNet(cellchat)

#保存文件
saveRDS(cellchat, file = "./data/temp/T-Duct_cellchat.all_tumor.rds")

##7. 单个信号通路####
cellchat <- readRDS(file = "./data/temp/T-Duct_cellchat_tumor.rds")
## Access all the signaling pathways showing significant communications
pathways.show.all <- cellchat@netP$pathways
pathways.show.all
# [1] "MHC-I"     "COLLAGEN"  "LAMININ"   "MIF"       "MK"        "CLEC"      "CD99"      "FN1"      
# [9] "GALECTIN"  "APP"       "SPP1"      "ADGRE5"    "MHC-II"    "CDH1"      "CDH"       "JAM"      
# [17] "VISFATIN"  "PARs"      "EGF"       "NECTIN"    "TGFb"      "CD46"      "LCK"       "TIGIT"    
# [25] "SEMA4"     "EPHA"      "CEACAM"    "DESMOSOME" "HSPG"      "GRN"       "EPHB"      "MPZ"      
# [33] "THBS"      "GDF"       "ALCAM"     "CD6"       "OCLN"      "AGRN"      "SEMA3"     "TNF"      
# [41] "CD96"      "CCL"       "PTN"       "PVR"       "SELPLG"    "NOTCH"     "ANGPTL"    "IFN-II"   
# [49] "ITGB2"     "SELL"      "CXCL"      "TRAIL"     "CD137"     "BMP"       "WNT"       "SEMA7"    
# [57] "VCAM"      "ICAM"      "GAS"       "SEMA6"     "PROS"      "VEGI"      "ESAM"  

# select one pathway
pathways.show <- c("TIGIT") #⭐T和Duct的交互
pathways.show <- c("EGF")#内部互相促进
pathways.show <- c("TNF")#⭐ OK
pathways.show <- c("SPP1")#
pathways.show <- c("LAMININ")#配对中有一个分别作用T和Duct，与tumor匹配
pathways.show <- c("MHC-II")#HLA-DMA - CD4,与tumor匹配
#舍弃

pathways.show <- c("CDH")
pathways.show <- c("MHC-I") #
pathways.show <- c("ADGRE5")#仅一个配对
pathways.show <- c("PARs")#无特色
pathways.show <- c("CD99")#仅一个配对
pathways.show <- c("APP")#仅一个配对
pathways.show <- c("COLLAGEN")
pathways.show <- c("MIF")
pathways.show <- c("FN1")
pathways.show <- c("TGFb")
pathways.show <- c("NOTCH")#
pathways.show <- c("WNT")#
pathways.show <- c("LAMININ")#
pathways.show <- c("EPHA")


###1. 环状图####
par(mfrow=c(1,1))
p1 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle",  color.use = NULL, sources.use = NULL, targets.use = NULL, idents.use = NULL)#label.edge= T,#显示贡献值

###2. 层次聚类图####
vertex.receiver = seq(1,9)
p2 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout ="hierarchy", vertex.receiver = vertex.receiver)
# 左半部分是自分泌相关信号，自己释放的信号作用于自己
# 右半部分就是展示的旁分泌信号
# 链接：https://www.jianshu.com/p/cf79a2b7703b

###3. 和弦图####
par(mfrow=c(1,1))
p3 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout ="chord")

par(mfrow=c(1,1))
group.cellType <- c(rep("CD4", 4), rep("CD8", 5),rep("Duct", 9), rep("NKT", 1) ) #grouping cell clusters into fibroblast, DC and TC cells
names(group.cellType) <- levels(cellchat@idents)
p4 <- netVisual_chord_cell(cellchat, signaling = pathways.show, group =group.cellType, title.name = paste0(pathways.show, " signaling network"))
p4
###4. 热图####
par(mfrow=c(1,1))
p5 <- netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
p5
###5. 柱状图#### 
# 配受体贡献列表
p6 <- netAnalysis_contribution(cellchat, signaling = pathways.show)
p6
#配受体贡献列表配受体贡献列表
pairLR <- extractEnrichedLR(cellchat, signaling = pathways.show, geneLR.return = FALSE)
pairLR

###6. 配受体作用情况####
LR.show <- pairLR[1,]
p71 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[2,]
p72 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[3,]
p73 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[4,]
p74 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[5,]
p75 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[6,]
p76 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[7,]
p77 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[8,]
p78 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[9,]
p79 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[10,]
p710 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[11,]
p711 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[12,]
p712 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[13,]
p713 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")


pdf("./data/output/T-Duct_cellchat_LAMININ_tumor.pdf",width =10 ,height = 10)
p1
p2
p3
p4
p5
p6
p71
p72
p73
p74
p75
p76
p77
p78
p79
p710
p711
p712
p713
dev.off()

###7. 配受体小提琴图####
pathways.show <- c("TIGIT")

#通路配受体表达强度，小提琴图
plotGeneExpression(cellchat, signaling = pathways.show, enriched.only =TRUE)

#计算并可视化网络中心性得分
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")# the slot 'netP' means the inferred intercellular communication network of signaling pathways

#使用热图将计算的中心性得分可视化，识别细胞群的主要信号作用
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

#在2D空间中可视化主要源和目标
netAnalysis_signalingRole_scatter(cellchat, signaling = NULL)

#识别对某些细胞组的传出或传入信号贡献最大的信号
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")
ht1 + ht2


##8. 气泡图-出发结束####
# sources.use, targets.use 对应序号 
levels(cellchat@idents)
# [1] "CD4_Tem"  "CD4_Th"   "CD4_Tn"   "CD4_Treg" "CD8_Tc"   "CD8_Te"   "CD8_Tem"  "CD8_Tex"  "CD8_Trm" 
# [10] "Group_1"  "Group_2"  "Group_3"  "Group_4"  "Group_5"  "Group_6"  "Group_7"  "Group_8"  "Group_9" 
# [19] "NKT" 

###1. 所有的信号通路####
# T出发
p1 <- netVisual_bubble(cellchat, sources.use = c(1:9,19), targets.use = c(10:18), remove.isolate = FALSE)#去掉空白行和空白列
# Duct出发
p2 <- netVisual_bubble(cellchat, sources.use = c(10:18), targets.use = c(1:9,19), remove.isolate = FALSE)
pdf("./data/output/T-Duct_cellchat_All_气泡图_tumor.pdf",width =15 ,height = 10)
p1
p2
dev.off()

###2. 看指定信号通路####
pairLR.use <- extractEnrichedLR(cellchat, signaling = c("TIGIT","EGF","TNF","SPP1","MIF","MHC-II"))

netVisual_bubble(cellchat, sources.use = c(1:9,19), targets.use = c(10:18), pairLR.use = pairLR.use, remove.isolate = TRUE)
netVisual_bubble(cellchat, sources.use = c(10:18), targets.use = c(1:9,19),  
                 # pairLR.use = pairLR.use, 
                 remove.isolate = TRUE)

###3. 看内部信号通路####
#CD4内部 CD4Treg起点
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(1:4), 
                 #signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), 
                 remove.isolate = FALSE)
#CD4内部 CD4Treg终点
netVisual_bubble(cellchat, targets.use = 4, sources.use = c(1:4), 
                 # signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), 
                 remove.isolate = FALSE)
#CD8内部 CD8Tex起点
netVisual_bubble(cellchat, sources.use = 8, targets.use = c(6:10), 
                 signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), remove.isolate = FALSE)
#CD8内部 CD8Tex终点
netVisual_bubble(cellchat, targets.use = 8, sources.use = c(6:10), 
                 signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), remove.isolate = FALSE)

#列出所有与细胞相关的配体受体⭐
netVisual_bubble(cellchat, targets.use = c("CD4_Treg","CD8_Tex","NKT"), pairLR.use = pairLR.use, 
                 remove.isolate = TRUE, sort.by.target = T)
netVisual_bubble(cellchat, sources.use =c("CD4_Treg","CD8_Tex","NKT"), pairLR.use = pairLR.use, remove.isolate = TRUE,sort.by.source = T)

netVisual_bubble(cellchat, 
                 sources.use = c("CD4_Treg","CD8_Tex","NKT"), 
                 targets.use = c("CD4_Treg","CD8_Tex","NKT"), 
                 pairLR.use = pairLR.use, 
                 remove.isolate = TRUE, 
                 sort.by.source = T, 
                 sort.by.target = T, 
                 sort.by.source.priority = F)#排序优先级
##9. 基因和信号通路####
###1. 和弦图####

netVisual_chord_gene(cellchat, sources.use = 4, targets.use = c(1:4), lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = 8, targets.use = c(5:10), lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = 10, targets.use = c(5:10), lab.cex = 0.5,legend.pos.y = 30)

netVisual_chord_gene(cellchat, sources.use = c(1:4), targets.use = 10, lab.cex =0.5,legend.pos.y = 40)
netVisual_chord_gene(cellchat, sources.use = c(5:10), targets.use = 8, lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = c(5:10), targets.use = 10, lab.cex = 0.5,legend.pos.y = 30)

# show all the significant interactions (L-R pairs) associated with certain signaling pathways
# netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("CCL","CXCL"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c("CD4_Treg","CD8_Tex","NKT"), targets.use = c("CD4_Treg","CD8_Tex","NKT"), signaling =  c("TIGIT","EGF","TNF","SPP1","LAMININ","MHC-II"),legend.pos.x = 8, lab.cex = 0.5)

# show all the significant signaling pathways from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), slot.name = "netP", legend.pos.x = 10)


##10. 通信模式####
#(A)识别和可视化分泌细胞的outgoing通信模式
# 推断图案的数量。
selectK(cellchat, pattern ="outgoing")
#当传出模式的数量为3时，Cophenetic和Silhouette值都开始突然下降。
nPatterns = 6 #选择合适的分组数
cellchat <- identifyCommunicationPatterns(cellchat, pattern ="outgoing", k = nPatterns)
p1 <- netAnalysis_river(cellchat, pattern ="outgoing")# river plot
p2 <- netAnalysis_dot(cellchat, pattern ="outgoing")# dot plot

#(B)识别和可视化目标细胞的incoming通信模式
selectK(cellchat, pattern = "incoming")
nPatterns = 6 #
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns)
p3 <- netAnalysis_river(cellchat, pattern = "incoming")# river plot
p4 <- netAnalysis_dot(cellchat, pattern = "incoming")# dot plot

pdf("./data/output/T-Duct_cellchat_All_通信模式_tumor.pdf",width =8 ,height = 7)
p1
p2
p3
p4
dev.off()


##11. 结构和功能相似性####
###1. 功能相似性Functional similarity####
cellchat <- computeNetSimilarity(cellchat, type = "functional")
# 报错
# Manifold learning of the signaling networks for a single dataset 
# Error in runUMAP(Similarity, min_dist = min_dist, n_neighbors = n_neighbors,  : 
#                    Cannot find UMAP, please install through pip (e.g. pip install umap-learn or reticulate::py_install(packages = 'umap-learn')).
# 解决办法
# library(reticulate)
# use_python('C:\\Users\\ZFB\\AppData\\Local\\Programs\\Python\\Python39\\',required = T)
# py_version() #看看是不是指定正确
# cellchat <- netEmbedding(cellchat, type = "functional") #未解决

#install through "pip install umap-learn"
#或者下面三條命令都可以運行
# cellchat <- netEmbedding(cellchat, umap.method='uwot',type ="functional")
# cellchat <- netEmbedding(cellchat, umap.method='umap-learn',type ="functional")
cellchat <- netClustering(cellchat, type = "functional",do.parallel=FALSE)#加上"do.parallel=FALSE"
netVisual_embedding(cellchat, type = "functional", label.size =3.5)
netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)

###2. 结构相似性Structure similarity####
cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")#或者下面這條命令
cellchat <- netEmbedding(cellchat, umap.method='uwot',type = "structural")
cellchat <- netClustering(cellchat, type = "structural",do.parallel=FALSE)#加上"do.parallel=FALSE"
netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
netVisual_embeddingZoomIn(cellchat, type = "structural", nCol = 2)

saveRDS(cellchat, file = "./data/temp/T-Duct_cellchat.all_tumor.rds")


#⭐ normal####
##1. 读入数据####
obj_all <- readRDS("./data/temp/T-Duct_cellchat.all_obj_all.rds")
obj_normal <- subset(obj_all, tech=="Normal")

##2. 转化为cellchat数据格式####
cellchat <- createCellChat(object = obj_normal, group.by = "labels", assay = "RNA")

##3. 设置参考数据库####
# 根据分析数据的物种，可选CellChatDB.human, 或者 CellChatDB.mouse 。通过showDatabaseCategory函数可以查看该数据库的情况
CellChatDB <- CellChatDB.human
showDatabaseCategory(CellChatDB)
# 展示参考数据库
dplyr::glimpse(CellChatDB$interaction)
# 使用默认所有的参考数据库
CellChatDB.use <- CellChatDB 
# 合并保存到cellchat内
cellchat@DB <- CellChatDB.use

# 对信号基因的表达数据取子集以节省计算成本
cellchat <- subsetData(cellchat)                         # 即使使用整个数据库，此步骤也是必要的
future::plan("multisession", workers = 10)                # 运行线程数
cellchat <- identifyOverExpressedGenes(cellchat)         #识别过表达基因
cellchat <- identifyOverExpressedInteractions(cellchat)  #识别过表达配体受体对

#细胞通信网络的推理
cellchat <- computeCommunProb(cellchat, type = "triMean") #triMean用于计算每个细胞组的平均基因表达。

# cellchat <- filterCommunication(cellchat, min.cells = 10)
# cellchat <- computeCommunProbPathway(cellchat)

##4. 推断细胞通讯网络####
#使用表达值推测细胞互作的概率，该步骤相对较耗时一些。
cellchat <- computeCommunProb(cellchat, raw.use = TRUE, population.size = TRUE) 
# Filter out the cell-cell communication if there are only few number of cells in certain cell groups
# 注1：raw.use = TRUE 表示使用raw数据，而不使用上一步projectData后的结果。
# 注2：在假设细胞数较多的群 往往比 细胞数较少的群发送更强的信号的前提下，当population.size = TRUE时候，CellChat可以在概率计算中考虑每个细胞群中细胞比例的影响。
cellchat <- filterCommunication(cellchat, min.cells = 10)



##5. 提取 保存结果####
#all the inferred cell-cell communications at the level of ligands/receptors
df.net <- subsetCommunication(cellchat)
write.csv(df.net, "./data/output/T-Duct_cellchat.all_normal.csv")

#获取显著的结果
df.net1 <- subsetCommunication(cellchat,slot.name = "netP")

#获取起始和结束的关系
levels(cellchat@idents)
df.net2 <- subsetCommunication(cellchat, sources.use = c("Epi"), targets.use = c("Fibroblast" ,"T")) 

#获取信号通路
df.net3 <- subsetCommunication(cellchat, signaling = c("CCL", "TGFb"))

##6. 计算cell-cell communication####
#计算每个信号通路相关的所有配体-受体相互作用的通信结果，结存存放在net 和 netP中 。
cellchat <- computeCommunProbPathway(cellchat)

#计算整合的细胞类型之间通信结果
cellchat <- aggregateNet(cellchat)

#保存文件
saveRDS(cellchat, file = "./data/temp/T-Duct_cellchat.all_normal.rds")

##7. 单个信号通路####
cellchat <- readRDS(file = "./data/temp/T-Duct_cellchat.all_normal.rds")
## Access all the signaling pathways showing significant communications
pathways.show.all <- cellchat@netP$pathways
pathways.show.all
# [1] "MHC-I"     "COLLAGEN"  "LAMININ"   "MIF"       "MK"        "CLEC"      "CD99"      "FN1"      
# [9] "GALECTIN"  "APP"       "SPP1"      "ADGRE5"    "MHC-II"    "CDH1"      "CDH"       "JAM"      
# [17] "VISFATIN"  "PARs"      "EGF"       "NECTIN"    "TGFb"      "CD46"      "LCK"       "TIGIT"    
# [25] "SEMA4"     "EPHA"      "CEACAM"    "DESMOSOME" "HSPG"      "GRN"       "EPHB"      "MPZ"      
# [33] "THBS"      "GDF"       "ALCAM"     "CD6"       "OCLN"      "AGRN"      "SEMA3"     "TNF"      
# [41] "CD96"      "CCL"       "PTN"       "PVR"       "SELPLG"    "NOTCH"     "ANGPTL"    "IFN-II"   
# [49] "ITGB2"     "SELL"      "CXCL"      "TRAIL"     "CD137"     "BMP"       "WNT"       "SEMA7"    
# [57] "VCAM"      "ICAM"      "GAS"       "SEMA6"     "PROS"      "VEGI"      "ESAM"  

# select one pathway
pathways.show <- c("TIGIT") #⭐T和Duct的交互
pathways.show <- c("EGF")#内部互相促进
pathways.show <- c("TNF")#⭐ OK
pathways.show <- c("SPP1")#
pathways.show <- c("LAMININ")#配对中有一个分别作用T和Duct，与tumor匹配
pathways.show <- c("MHC-II")#HLA-DMA - CD4,与tumor匹配
#舍弃

pathways.show <- c("CDH")
pathways.show <- c("MHC-I") #
pathways.show <- c("ADGRE5")#仅一个配对
pathways.show <- c("PARs")#无特色
pathways.show <- c("CD99")#仅一个配对
pathways.show <- c("APP")#仅一个配对
pathways.show <- c("COLLAGEN")
pathways.show <- c("MIF")
pathways.show <- c("FN1")
pathways.show <- c("TGFb")
pathways.show <- c("NOTCH")#
pathways.show <- c("WNT")#
pathways.show <- c("LAMININ")#
pathways.show <- c("EPHA")


###1. 环状图####
par(mfrow=c(1,1))
p1 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle",  color.use = NULL, sources.use = NULL, targets.use = NULL, idents.use = NULL)#label.edge= T,#显示贡献值

###2. 层次聚类图####
vertex.receiver = seq(1,9)
p2 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout ="hierarchy", vertex.receiver = vertex.receiver)
# 左半部分是自分泌相关信号，自己释放的信号作用于自己
# 右半部分就是展示的旁分泌信号
# 链接：https://www.jianshu.com/p/cf79a2b7703b

###3. 和弦图####
par(mfrow=c(1,1))
p3 <- netVisual_aggregate(cellchat, signaling = pathways.show, layout ="chord")

par(mfrow=c(1,1))
group.cellType <- c(rep("CD4", 4), rep("CD8", 5),rep("Duct", 9), rep("NKT", 1) ) #grouping cell clusters into fibroblast, DC and TC cells
names(group.cellType) <- levels(cellchat@idents)
p4 <- netVisual_chord_cell(cellchat, signaling = pathways.show, group =group.cellType, title.name = paste0(pathways.show, " signaling network"))
p4
###4. 热图####
par(mfrow=c(1,1))
p5 <- netVisual_heatmap(cellchat, signaling = pathways.show, color.heatmap = "Reds")
p5
###5. 柱状图#### 
# 配受体贡献列表
p6 <- netAnalysis_contribution(cellchat, signaling = pathways.show)
p6
#配受体贡献列表配受体贡献列表
pairLR <- extractEnrichedLR(cellchat, signaling = pathways.show, geneLR.return = FALSE)
pairLR

###6. 配受体作用情况####
LR.show <- pairLR[1,]
p71 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[2,]
p72 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[3,]
p73 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[4,]
p74 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[5,]
p75 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[6,]
p76 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[7,]
p77 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[8,]
p78 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[9,]
p79 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[10,]
p710 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[11,]
p711 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[12,]
p712 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")
LR.show <- pairLR[13,]
p713 <- netVisual_individual(cellchat, signaling = pathways.show, pairLR.use = LR.show, layout = "circle")


pdf("./data/output/T-Duct_cellchat_MHC-II_normal.pdf",width =10 ,height = 10)
p1
p2
p3
p4
p5
p6
p71
p72
p73
p74
p75
p76
p77
p78
p79
p710
p711
p712
p713
dev.off()

###7. 配受体小提琴图####
pathways.show <- c("TIGIT")

#通路配受体表达强度，小提琴图
plotGeneExpression(cellchat, signaling = pathways.show, enriched.only =TRUE)

#计算并可视化网络中心性得分
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")# the slot 'netP' means the inferred intercellular communication network of signaling pathways

#使用热图将计算的中心性得分可视化，识别细胞群的主要信号作用
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show, width = 8, height = 2.5, font.size = 10)

#在2D空间中可视化主要源和目标
netAnalysis_signalingRole_scatter(cellchat, signaling = NULL)

#识别对某些细胞组的传出或传入信号贡献最大的信号
ht1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing")
ht2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")
ht1 + ht2


##8. 气泡图-出发结束####
# sources.use, targets.use 对应序号 
levels(cellchat@idents)
# [1] "CD4_Tem"  "CD4_Th"   "CD4_Tn"   "CD4_Treg" "CD8_Tc"   "CD8_Te"   "CD8_Tem"  "CD8_Tex"  "CD8_Trm" 
# [10] "Group_1"  "Group_2"  "Group_3"  "Group_4"  "Group_5"  "Group_6"  "Group_7"  "Group_8"  "Group_9" 
# [19] "NKT" 

###1. 所有的信号通路####
# T出发
p1 <- netVisual_bubble(cellchat, sources.use = c(1:9,19), targets.use = c(10:18), remove.isolate = FALSE)#去掉空白行和空白列
# Duct出发
p2 <- netVisual_bubble(cellchat, sources.use = c(10:18), targets.use = c(1:9,19), remove.isolate = FALSE)
pdf("./data/output/T-Duct_cellchat_All_气泡图_tumor.pdf",width =15 ,height = 10)
p1
p2
dev.off()

###2. 看指定信号通路####
pairLR.use <- extractEnrichedLR(cellchat, signaling = c("TIGIT","EGF","TNF","SPP1","MIF","MHC-II"))

netVisual_bubble(cellchat, sources.use = c(1:9,19), targets.use = c(10:18), pairLR.use = pairLR.use, remove.isolate = TRUE)
netVisual_bubble(cellchat, sources.use = c(10:18), targets.use = c(1:9,19),  
                 # pairLR.use = pairLR.use, 
                 remove.isolate = TRUE)

###3. 看内部信号通路####
#CD4内部 CD4Treg起点
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(1:4), 
                 #signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), 
                 remove.isolate = FALSE)
#CD4内部 CD4Treg终点
netVisual_bubble(cellchat, targets.use = 4, sources.use = c(1:4), 
                 # signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), 
                 remove.isolate = FALSE)
#CD8内部 CD8Tex起点
netVisual_bubble(cellchat, sources.use = 8, targets.use = c(6:10), 
                 signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), remove.isolate = FALSE)
#CD8内部 CD8Tex终点
netVisual_bubble(cellchat, targets.use = 8, sources.use = c(6:10), 
                 signaling = c("MHC-I","MIF","CLEC","COLLAGEN","LCK","ITGB2"), remove.isolate = FALSE)

#列出所有与细胞相关的配体受体⭐
netVisual_bubble(cellchat, targets.use = c("CD4_Treg","CD8_Tex","NKT"), pairLR.use = pairLR.use, 
                 remove.isolate = TRUE, sort.by.target = T)
netVisual_bubble(cellchat, sources.use =c("CD4_Treg","CD8_Tex","NKT"), pairLR.use = pairLR.use, remove.isolate = TRUE,sort.by.source = T)

netVisual_bubble(cellchat, 
                 sources.use = c("CD4_Treg","CD8_Tex","NKT"), 
                 targets.use = c("CD4_Treg","CD8_Tex","NKT"), 
                 pairLR.use = pairLR.use, 
                 remove.isolate = TRUE, 
                 sort.by.source = T, 
                 sort.by.target = T, 
                 sort.by.source.priority = F)#排序优先级
##9. 基因和信号通路####
###1. 和弦图####

netVisual_chord_gene(cellchat, sources.use = 4, targets.use = c(1:4), lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = 8, targets.use = c(5:10), lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = 10, targets.use = c(5:10), lab.cex = 0.5,legend.pos.y = 30)

netVisual_chord_gene(cellchat, sources.use = c(1:4), targets.use = 10, lab.cex =0.5,legend.pos.y = 40)
netVisual_chord_gene(cellchat, sources.use = c(5:10), targets.use = 8, lab.cex = 0.5,legend.pos.y = 30)
netVisual_chord_gene(cellchat, sources.use = c(5:10), targets.use = 10, lab.cex = 0.5,legend.pos.y = 30)

# show all the significant interactions (L-R pairs) associated with certain signaling pathways
# netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), signaling = c("CCL","CXCL"),legend.pos.x = 8)
netVisual_chord_gene(cellchat, sources.use = c("CD4_Treg","CD8_Tex","NKT"), targets.use = c("CD4_Treg","CD8_Tex","NKT"), signaling =  c("TIGIT","EGF","TNF","SPP1","LAMININ","MHC-II"),legend.pos.x = 8, lab.cex = 0.5)

# show all the significant signaling pathways from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
netVisual_chord_gene(cellchat, sources.use = c(1,2,3,4), targets.use = c(5:11), slot.name = "netP", legend.pos.x = 10)


##10. 通信模式####
#(A)识别和可视化分泌细胞的outgoing通信模式
# 推断图案的数量。
selectK(cellchat, pattern ="outgoing")
#当传出模式的数量为3时，Cophenetic和Silhouette值都开始突然下降。
nPatterns = 6 #选择合适的分组数
cellchat <- identifyCommunicationPatterns(cellchat, pattern ="outgoing", k = nPatterns)
p1 <- netAnalysis_river(cellchat, pattern ="outgoing")# river plot
p2 <- netAnalysis_dot(cellchat, pattern ="outgoing")# dot plot

#(B)识别和可视化目标细胞的incoming通信模式
selectK(cellchat, pattern = "incoming")
nPatterns = 6 #
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = nPatterns)
p3 <- netAnalysis_river(cellchat, pattern = "incoming")# river plot
p4 <- netAnalysis_dot(cellchat, pattern = "incoming")# dot plot

pdf("./data/output/T-Duct_cellchat_All_通信模式_tumor.pdf",width =8 ,height = 7)
p1
p2
p3
p4
dev.off()


##11. 结构和功能相似性####
###1. 功能相似性Functional similarity####
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
#install through "pip install umap-learn"
#或者下面三條命令都可以運行
cellchat <- netEmbedding(cellchat, umap.method='uwot',type ="functional")
cellchat <- netEmbedding(cellchat, umap.method='umap-learn',type ="functional")
cellchat <- netClustering(cellchat, type = "functional",do.parallel=FALSE)#加上"do.parallel=FALSE"
netVisual_embedding(cellchat, type = "functional", label.size =3.5)
netVisual_embeddingZoomIn(cellchat, type = "functional", nCol = 2)

###2. 结构相似性Structure similarity####
cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")#或者下面這條命令
cellchat <- netEmbedding(cellchat, umap.method='uwot',type = "structural")
cellchat <- netClustering(cellchat, type = "structural",do.parallel=FALSE)#加上"do.parallel=FALSE"
netVisual_embedding(cellchat, type = "structural", label.size = 3.5)
netVisual_embeddingZoomIn(cellchat, type = "structural", nCol = 2)

saveRDS(cellchat, file = "./data/temp/T-Duct_cellchat.all_tumor.rds")

saveRDS(cellchat, file = "./data/temp/T-Duct_cellchat.all_normal.rds")




#未使用 Part 2. ####
#Comparative analysis of cell-cell communication from pairs of scRNA-seq datasets
cellchat.LS <- readRDS("cellchat_humanSkin_LS.rds")
cellchat.NL <- readRDS("cellchat_humanSkin_NL.rds")
cellchat.NL <- updateCellChat(cellchat)
cellchat.LS <- updateCellChat(cellchat.LS)
object.list <- list(NL = cellchat.NL, LS = cellchat.LS)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))
save(object.list, file = "cellchat_object.list_humanSkin_NL_LS.RData")
save(cellchat, file = "cellchat_merged_humanSkin_NL_LS.RData")
gg1 <- compareInteractions(cellchat, show.legend = F, group = c (1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c (1,2), measure = "weight")
gg1 + gg2

#(A) Circle plot showing differential number of interactions or interaction strength among different cell populations across two datasets
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T)#install igraph  version1.3.5  #https://blog.csdn.net/m0_55681975/article/details/131965802
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")

#(B) Heatmap showing differential number of interactions or interaction strength among different cell populations across two datasets
gg1 <- netVisual_heatmap(cellchat)
gg2 <- netVisual_heatmap(cellchat, measure = "weight")
gg1 + gg2

#(C) Circle plot showing the number of interactions or interaction strength among different cell populations across multiple datasets
weight.max <- getMaxWeight(object.list, attribute = c("idents","count"))
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$count, weight.scale = T,
                   label.edge= F, edge.weight.max = weight.max[2], edge.width.max =
                     12, title.name = paste0("Number of interactions - ", names(object.list)[i]))}

#(D) Circle plot showing the differential number of interactions or interaction strength among coarse cell types
# Here, CellChat categorize the cell populations into three cell types, and then re-merge the list of CellChat objects.
group.cellType <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4))
group.cellType <- factor(group.cellType, levels = c("FIB", "DC","TC"))
object.list <- lapply(object.list, function(x) {mergeInteractions (x, group.cellType)})
cellchat <- mergeCellChat(object.list, add.names = names(object.list))

weight.max <- getMaxWeight(object.list, slot.name = c("idents", "net", "net"), attribute = c("idents","count", "count.merged"))
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$count.merged, weight.scale = T,
                   label.edge= T, edge.weight.max = weight.max[3], edge.width.max = 12, 
                   title.name = paste0("Number of interactions - ", names(object.list)[i]))}

# Similarly, CellChat can also show the differential number of interactions or interaction strength between any two cell types using circle plot.
par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "count.merged", label.edge = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight.merged", label.edge = T)

#(A) Identify cell populations with significant changes in sending or receiving signals
num.link <- sapply(object.list, function(x) {rowSums(x@net$count) + colSums(x@net$count)-diag(x@net$count)})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets
gg <- list()
for (i in 1:length(object.list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]],
                                               title = names(object.list)[i], weight.MinMax = weight.MinMax)}
patchwork::wrap_plots(plots = gg)

#(B) Identify the signaling changes of specific cell populations
gg1 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use  = "Inflam. DC", signaling.exclude = "MIF")
gg2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "cDC1", signaling.exclude = c("MIF"))
patchwork::wrap_plots(plots = list(gg1,gg2))

#remotes::install_github("hafen/rminiconda")
#py <- rminiconda::find_miniconda_python("my_python")
#reticulate::use_python("/home/yueli/.local/share/r-miniconda/envs/r-reticulate/bin/python", required = TRUE)
#reticulate::py_install(packages = 'umap-learn')

cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional") 
#cellchat <- netEmbedding(cellchat, umap.method="uwot",type = "functional") 
#cellchat <- netEmbedding(cellchat, umap.method="umap-learn",type = "functional")

cellchat <- netClustering(cellchat, type = "functional",do.parallel=FALSE)#
#cellchat <- netClustering(cellchat, method="umap-learn", type = "functional")#not work
#cellchat <- netClustering(cellchat, method="uwot", type = "functional")#not work
netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)
netVisual_embeddingPairwiseZoomIn(cellchat, type = "functional", nCol = 2)


rankSimilarity(cellchat, type = "functional")
gg1 <- rankNet(cellchat, mode = "comparison", stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", stacked = F, do.stat = TRUE)
gg1 + gg2

#(B) Compare outgoing (or incoming) signaling patterns associated with each cell population
i = 1
# combining all the identified signaling pathways from different datasets
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "outgoing", 
                                        signaling = pathway.union, title = names(object.list)[i], width = 5, height= 6)
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "outgoing", 
                                        signaling = pathway.union, title = names(object.list)[i+1], width = 5, height = 6)
draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

# Compare the communication probabilities from certain cell groups to other cell groups
netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11),comparison= c(1, 2), angle.x = 45)
# Identify the up-regulated (increased) and down-regulated (decreased) signaling ligand-receptor pairs in one dataset compared to the other dataset.                 
gg1 <- netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11), comparison = c(1, 2), max.dataset = 2, title.name = "Increased signaling in LS",
                        angle.x = 45, remove.isolate = T)
gg2 <- netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11), comparison = c(1, 2), max.dataset = 1, title.name = "Decreased signaling in LS",
                        angle.x = 45, remove.isolate = T)
gg1 + gg2

# define a positive dataset, i.e., the dataset with positive fold change against the other dataset
pos.dataset = "LS"
# define a char name used for storing the results of differential expression analysis
features.name = pos.dataset
# perform differential expression analysis
cellchat <- identifyOverExpressedGenes(cellchat, group.dataset = "datasets", pos.dataset = pos.dataset, 
                                       features.name = features.name, only.pos = FALSE, thresh.pc = 0.1, thresh.fc = 0.1, thresh.p = 1)
# map the results of differential expression analysis onto the inferred cell-cell communications to easily manage/subset the ligand-receptor pairs of interest
net <- netMappingDEG(cellchat, features.name = features.name)
# extract the ligand-receptor pairs with upregulated ligands in LS
net.up <- subsetCommunication(cellchat, net = net, datasets = "LS",ligand.logFC = 0.2, receptor.logFC = NULL)
# extract the ligand-receptor pairs with upregulated ligands and upregulated recetptors in NL, i.e.,downregulated in LS
net.down <- subsetCommunication(cellchat, net = net, datasets = "NL",ligand.logFC = -0.1, receptor.logFC = -0.1)
# do further deconvolution to obtain the individual signaling genes
gene.up <- extractGeneSubsetFromPair(net.up, cellchat)
gene.down <- extractGeneSubsetFromPair(net.down, cellchat)
# Users can also find all the significant outgoing/incoming/both signaling according to the customized DEG features and cell groups of interest
#df <- findEnrichedSignaling(object.list[[2]], features = c("CCL19", "CXCL12"), idents = c("Inflam. FIB", "COL11A1+ FIB"), pattern ="outgoing")  #It is availabe in CellChat version 2.1.1 
#Error in findEnrichedSignaling(object.list[[2]], features = c("CCL19",  : 
#                                                               could not find function "findEnrichedSignaling"pairLR.use.up = net.up[, "interaction_name", drop = F]

#(A) Bubble plot
pairLR.use.up = net.up[, "interaction_name", drop = F]
gg1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = 4, targets.use = c(5:11), comparison = c(1, 2), angle.x = 90, remove.isolate =  T,title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
pairLR.use.down = net.down[, "interaction_name", drop = F]
gg2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = 4, targets.use = c(5:11), comparison = c(1, 2), angle.x = 90, remove.isolate = T,title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))
gg1 + gg2

#Chord diagram
par(mfrow = c(1,2), xpd=TRUE)
netVisual_chord_gene(object.list[[2]], sources.use = 4, targets.use = c(5:11),
                     slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name =
                       paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 4, targets.use = c(5:11),
                     slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name
                     = paste0("Down-regulated signaling in ", names(object.list)[2]))

#Wordcloud plot
# visualize the enriched ligands in the first condition
computeEnrichmentScore(net.down, species = 'human')
# visualize the enriched ligands in the second condition
computeEnrichmentScore(net.up, species = 'human')

#(A) Circle plot
pathways.show <- c("CXCL")
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name =
                        paste(pathways.show, names(object.list)[i]))}
#Heatmap plot
pathways.show <- c("CXCL")
par(mfrow = c(1,2), xpd=TRUE)
ht <- list()
for (i in 1:length(object.list)) {
  ht[[i]] <- netVisual_heatmap(object.list[[i]], signaling = pathways.show, color.heatmap = "Reds",title.name = paste(pathways.show, "signaling ",names(object.list)[i]))
}
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))

#Visualize gene expression distribution.
cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level
plotGeneExpression(cellchat, signaling = "CXCL", split.by = "datasets", colors.ggplot = T)

save(object.list, file = "cellchat_object.list_humanSkin_NL_LS.RData")
save(cellchat, file = "cellchat_merged_humanSkin_NL_LS.RData")

#https://ndownloader.figshare.com/files/25957094
#https://ndownloader.figshare.com/files/25957634

cellchat.E13 <- readRDS("cellchat_embryonic_E13.rds")
cellchat.E13 <- updateCellChat(cellchat.E13)
cellchat.E14 <- readRDS("cellchat_embryonic_E14.rds")
cellchat.E14 <- updateCellChat(cellchat.E14)

group.new = levels(cellchat.E14@idents) # Define the cell labels to lift up
cellchat.E13 <- liftCellChat(cellchat.E13, group.new)
object.list <- list(E13 = cellchat.E13, E14 = cellchat.E14)
cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)

# Circle plot
pathways.show <- c("WNT")
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[i]))}

save(object.list, file = "cellchat_object.list_embryonic_E13_E14.RData")
save(cellchat, file = "cellchat_merged_embryonic_E13_E14.RData")

#https://figshare.com/articles/dataset/10X_visium_data_for_spatial-informed_cell-cell_communication/23621151
load("visium_mouse_cortex_annotated.RData")

# Gene expression data
data.input = GetAssayData(visium.brain, slot = "data", assay = "SCT") # normalized data matrix
# User assigned cell labels
meta = data.frame(labels = Idents(visium.brain), row.names = names(Idents(visium.brain)))
# Spatial locations of spots from full (NOT high/low) resolution images
spatial.locs = GetTissueCoordinates(visium.brain, scale = NULL, cols = c("imagerow", "imagecol"))
# Scale factors and spot diameters of the full resolution images

#https://figshare.com/articles/dataset/spatial_imaging_data_for_the_10X_visium_brain_dataset/23709726
scale.factors = jsonlite::fromJSON('scalefactors_json.json')
scale.factors = list(spot.diameter = 65, spot = scale.factors$spot_diameter_fullres)

#Create a CellChat object
cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels", datatype = "spatial", coordinates = spatial.locs, scale.factors = scale.factors)
CellChatDB <- CellChatDB.mouse
# use a subset of CellChatDB for cell-cell communication analysis
CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling") # use Secreted Signaling
#CellChatDB.use <- CellChatDB # simply use the default CellChatDB to use all CellChatDB for cell-cell communication analysis
# set the used database in the object
cellchat@DB <- CellChatDB.use

#Identify over-expressed ligands or receptors.
cellchat <- subsetData(cellchat) # This step is necessary even ifusing the whole database
future::plan("multisession", workers = 4) # do parallel
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

#cellchat <- computeCommunProb(cellchat, type = "truncatedMean", trim = 0.1, distance.use = TRUE, interaction.range = 250, scale.distance = 0.01)
#long long time
cellchat <- computeCommunProb(cellchat, type = "truncatedMean", trim = 0.1, distance.use = TRUE,  scale.distance = 0.01)#long long time
# Filter the cell-cell communication
cellchat <- filterCommunication(cellchat, min.cells = 10)
# Infer the cell-cell communication at a signaling pathway level
cellchat <- computeCommunProbPathway(cellchat)
# Calculate the aggregated cell-cell communication network
cellchat <- aggregateNet(cellchat)
saveRDS(cellchat, file = "cellchat_visium_mouse_cortex.Rds")

groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1,2), xpd=TRUE)
netVisual_circle(cellchat@net$count, vertex.weight = rowSums(cellchat@net$count), weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(cellchat@net$weight, vertex.weight = rowSums(cellchat@net$weight), weight.scale = T, label.edge= F, title.name ="Interaction weights/strength")
pathways.show <- c("CXCL")
# Circle plot
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout = "circle")
# Spatial plot
par(mfrow=c(1,1))
netVisual_aggregate(cellchat, signaling = pathways.show, layout= "spatial", edge.width.max = 2, vertex.size.max = 1, alpha.image = 0.2, vertex.label.cex = 3.5)

#Updating the ligand-receptor interaction database CellChatDB

#https://htmlpreview.github.io/?https://github.com/jinworks/CellChat/blob/master/tutorial/Update-CellChat